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We present techniques for long-term, stable, and accurate evolutions of multiple-black-hole space- 
times using the 'moving puncture' approach with fourth- and eighth-order finite difference stencils. 
We use these techniques to explore configurations of three black holes in a hierarchical system con- 
sisting of a third black hole approaching a quasi-circular black-hole binary, and find that, depending 
on the size of the binary, the resulting encounter may lead to a prompt merger of all three black holes, 
production of a highly elliptical binary (with the third black hole remaining unbound), or disruption 
of the binary (leading to three free black holes). We also analyze the classical Burrau three-body 
problem using full numerical evolutions. In both cases, we find behaviors distinctly difi'erent from 
Newtonian predictions, which has important implications for N-body black-hole simulations. For 
our simulations we use approximate analytic initial data. We find that the eighth-order stencils sig- 
nificantly reduce the numerical errors for our choice of grid sizes, and that the approximate initial 
data produces the expected waveforms for black-hole binaries with modest initial separations. 
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I. INTRODUCTION 



The recent dramatic breakthroughs in the numerical 
techniques to evolve black-hole-binary spacetimes [EH, Hi 
has led to rapid advancements in our understanding of 
black-hole physics. Notable among these advancements 
are developments in mathematical relativity, including 
systems of PDEs and gauge choices d, the explo- 
ration of the validity of the cosmic censorship conjcc- 
ture iQ, |7|1 1 and the application of isolated horizon formu- 
lae [sllglliol. fill m, . There are man y e xciting new 
results on recoil velocities [l|Jl3JllJl^^ 



M M M M is m, [I^niiillialHinsI^ 

post-Newtonian (PN) and numerical waveform compar- 
iso ns fss]. [sgI [stI. 38 |, modeling of the remnant spin [gI, 
[Tol . Illl.l29 | . l39ll40l. 41 1, new studies of eccentric black hole 
binaries 7l.l42|.|43l a nd pro ducing waveforms for matched 
filtering [3, [H, [M [13, lH, ST 



I. I50l|. In particular, the re- 
cent discovery of very large merger recoil kicks for black- 
hole binaries with spins in the orbital plane, which was 
ori gina lly inferred from the results in ^21i] , then observed 
in [24| . and determined to have a maximum value of 
4000 kms~^ in '2^, has had a great impact in the as- 
trophysical community, with several groups now seeking 
for observational traces of such high speed holes as the 
byproduct of galaxy collisions H, |5l| . 



Three-body and four-body interactions are expected 
to be common in globular clusters [H, [s^ and in galac- 
tic cores hosting supermassive black holes (when stellar- 
mass-black-hole-binary systems interact with the super- 
massive black hole). Hierarchical triplets of supermassive 
black holes might also be formed in galactic nuclei under- 
going sequential mergers |,54i , i55i] . The gravitational wave 
emission from such systems was recently estimated using 
post-Newtonian techniques [56|. The recent discovery of 



a probable triple quasar |57j] (with estimated masses of 
50 X lO^M©, 100 X IO^Mq and 500 x lO^Mg and pro- 
jected separations of between 30 and 50 kpcs) indicates 
that hierarchical supermassive-three-black-hole systems 
are possible. Triple stars and black holes are much more 
common in globular clusters [ssj . and galactic disks. The 
closest star to the solar system, Alpha Centauri, is a 
triple system, as is Polaris and HD 188753. 

On the theoretical side, even in the Newtonian theory 
of gravity, the three-body problem is much more compli- 
cated than the two body one, and is generically chaotic. 
In Ref. we established that our numerical formalism 
is able to handle the evolution of three fully rclativistic 
bodies. In this paper we continue our quest to expand 
our understanding of the multiple-black-hole problem by 
studying a sequences of black-hole-binary — third-black- 
hole configurations, with the third black hole intersect- 
ing the binary along the axis of rotation, to explore the 
influence of a third black hole on the binary dynamics. 
We also examine the classical Burrau three-body problem 
and find that the dynamics can differ dramatically from 
the Newtonian prediction, with important consequences 
for N-body black-hole simulations. 

This paper is organized as follows: in Sec. |TT] we de- 
scribe the techniques used for multiple black-hole evolu- 
tions, including eighth-order techniques and approximate 
initial data. In Sec. IIIII we show the results from multi- 
ple black-hole evolutions of three families of hierarchical 
configurations consisting of a third black hole interacting 
with a black-hole binary along the binary's axis of rota- 
tion, as well as the classical Burrau three-body problem. 
We conclude our analysis in Sec. IIVI where we discuss 
the implications of our numerical results to black-hole 
astrophysics and N-body simulations. We provide the 
techniques to generate 2PN quasi-circular orbits of our 
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hierarchical three-black- hole systems in Appendix [Xl 



where 



II. TECHNIQUES 

We evolve the three-black-hole-data-sets using the 
LazEv [59| implementation of the 'moving puncture ap- 
proach' [2j_|3|. In our version of the moving puncture 
approach 0| we replace the BSSN [g^, [6l|, |62| confor- 
mal exponent which has logarithmic singularities at 
the punctures, with the initially field x — exp(— 4(/)). 
This new variable, along with the other BSSN variables, 
will remain finite provided that one uses a suitable choice 
for the gauge. We use the Carpet driver to provide a 
'moving boxes' style mesh refinement. In this approach 
refined grids of fixed size are arranged about the coordi- 
nate centers of each hole. The Carpet code then moves 
these fine grids about the computational domain by fol- 
lowing the trajectories of the three black holes. We use 
AHFinderDirect [ill to locate apparent horizons. 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with 
a modified l-|-log lapse and a modified Gamma-driver 
shift condition [2, l65l |. and an initial lapse a{t = 0) = 
2/(1-1- V'sl) (where tpBL is the Brill-Lindquist confor- 
mal factor discussed below). The lapse and shift are 
evolved with {pt - P^di)a = -2aK, dtfi" = S"^, and 
dtB"^ = i/AdtT'^ — rjB°^ . These gauge conditions require 
careful treatment of Xi the inverse of the three-metric 
conformal factor, near the pun cture in order for the sys- 
tem to remain stable gSi^. As was shown in Ref . Ijl , 
this choice of gauge leads to a strongly hyperbolic evolu- 
tion system provided that the shift does not become too 
large. Unless otherwise noted, we use a standard choice 
oi rj — &/AI for the simulations below. 



A. Approximate Initial Data 

We use the puncture approach [gS^ to compute initial 
data (See also [H, for similar approaches). In this 
approach the 3-metric on the initial slice has the form 
lab = [i'BL + u)^<5ah, whcrc ipBL is the Brill-Lindquist 
conformal factor, 5ab is the Euclidean metric, and u is 
(at least) on the punctures. The Brill-Lindquist con- 
formal factor is given by — '^i/(27'i)i where 
n is the total number of 'punctures', is the mass pa- 
rameter of puncture i {nii is not the horizon mass associ- 
ated with puncture i), and is the coordinate distance 
to puncture i. We solve for u using the approximate so- 
lutions given in Refs. [zEliilii, 

l74j . with the addition 
of a cross-term given below. We will use the notation 
of [7l| to write these approximate solutions. For a spin- 
ningpuncture, u has the form (Eqs. (4. 23), (4. 26), (4. 27) 
of 123): 
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(2a) 
(2b) 



J' = AJ/rn^, m is the mass parameter, I = 1/(1 + TZ), 
n = 2r/m,iJ,j = J-f, andP2(a;) = (3a:^-l)/2is the Leg- 
endre polynomial of degree 2. For a boosted puncture, u 
has the form (Eqs. (4.44), (4.48), (4.49) of 71]): 



where 
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21^ -f 21^ - ^4 + 
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(3) 



(4a) 



mHul = 15^ + 132^2 + 53^^ + 96^* + 82f^ + 
84fV7^ + 841n(^)/7^^ 



(4b) 



V = 2P/m and fip = P ■ f. If the puncture is both 
boosted and spinning then there is a cross-term 



Uc = iv X J)- 7^(l + 571 + lonyyso 



(5) 



For a spinning boosted puncture the approximate solu- 
tion has the form 



Up + uj - 



+ 0{V^) + 0{J^) 
0{V^J^) + 0{V^J), 



(6) 



where the error terms linear in V and J7 only occur when 
V X J ^ We obtain solutions for multiple punctures 
by superposition. The error in this approximation scales 
as the inverse of the distance squared between punctures. 

In this approximation the total ADM mass, linear mo- 
mentum, and angular momentum are given by: 



Madm = ^m, -f -j2/m3 + -p2/^^^ (7) 

i 

-PaDM = ^-Pij (8) 



•^ADM — y~^('^i 



xP.)- 



(9) 



UJ{^,^iJ) = {4 + uiTz^P2{fij)) , 



(1) 



We test the reliability of this approximate initial data 
by solving for equal-mass, non-spinning, quasi-circular 
binaries using both the TwoPunctures thorn 
(which is restricted to problems involving two punctures) 
and Eq. We use the 3PN equations of motion to ob- 
tain position and momentum parameters for two black 
holes in a binary with unit mass and a specified orbital 
frequency. We then choose puncture mass parameters 
(the same for both holes) such that the total ADM mass, 
as calculated using the TwoPunctures thorn, is IM. 
For our test case we choose a binary with orbital fre- 
quency Mft = 0.035 (which performs at-least three or- 
bits prior to merger) and evolve using a relatively coarse 
resolution of /i = M/32 (on the finest grid). As shown 
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FIG. 1: A comparison of the {£ — 2,m — 2) modes of 1/14 
at r = 40M for an equal-mass, non-spinning, quasi-circular 
binary of mass IM and orbital frequency Mfi = 0.035, pro- 
duced using the TwoPunctures thorn (labeled 'Exact') and 
the approximate data (labeled 'Approx'). In both cases we 
used identical values for the initial data parameters (puncture 
mass, puncture locations, puncture momenta). The approxi- 
mate data has been translated by St = 49. 5M and multiplied 
by a constant phase factor of exp(— 3.8289 i). The third plot 
('Approx_Renorm') shows the waveform using the approxi- 
mate data with the puncture masses rescaled such that the 
ADM mass as given by Eq. is IM. This latter plot has 
been translated by St = — 62A/ and multiplied by a phase 
factor of exp(4.5772 j). The inset shows an expanded view 
of the orbital motion. Note the excellent agreement in the 
waveforms. 
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B. Eighth-Order Finite Differencing 



We use both a new eighth-order spatial finite differ- 
encing algorithm and the standard fourth-order finite 
differencing used in our previous papers. Our eighth- 
order scheme extends the sixth-order scheme described in 
Ref. [t^]. As in [t^, we use a fourth-order Runge-Kutta 
time integrator and a second-order in time, fifth-order 
in space prolongation operator. Centered first spatial 
derivatives have the form: 



- 3/,+4 + 32/,+3 - 168/,+2 + 672/,+i) 
/{SAOdx) (10) 



(we suppress the other two indices in Eqs. pO)l - (fT3|) ). 
while for advection derivatives we adjust the center of 
the stencil by one point. The downward pointing stencil 
has the form: 



5x/, = (-3/,_5 + 30/,_4-140/,_3 + 420/,_2 - 1050/,;_i 
+378/, + 420/,+! - 60/,+2 + S/.+a) 
/{SAOdx), (11) 



while the upward pointing stencil has the form: 



in Fig. [U we find that the two waveforms agree after 
the first half-cycle of orbital motion (wc corrected for a 
time translation and phase difference in the waveform us- 
ing the techniques of [H, [ill). The phase difference and 
time translation between the waveforms produced by the 
two methods seems to be a result of the normalization of 
the approximate initial data (i.e. the puncture mass pa- 
rameters). When we normalize the approximate initial 
data such that the total approximate ADM mass (i.e. 
Eq. d?])) is IM (it is only exactly IM for the TwoPunc- 
tures data) , then the binary takes longer to inspiral and 
the phase difference and time translation are in the oppo- 
site direction. We find, based on a linear interpolation of 
the phase difference and time translation versus puncture 
mass, that there is a puncture mass choice where both 
the phase difference and time translation are zero. The 
puncture mass that gives an ADM mass of IM for the 
exact initial data is nip = 0.4891M, the puncture mass 
that gives and ADM mass of IM for the approximate 
data is nip = 0.4846M, and the interpolated puncture 
mass that gives phase (and time) agreement between the 
exact and approximate waveforms is nip = 0.4871M. 



- (+3/,+5 - 30/,+4 + 140/,+3 - 420/,+2 + 1050/,+i 
-378/, - 420/,_i -f- 60/,_2 - 5/,_3) 
/(840da;). (12) 



We use standard centered differencing for the second- 
spatial derivatives. The dxx, dyy, and dzz derivatives 
have the form: 



( - 9/,_4 + 128/,_3 - 1008/,_2 + 8064/i_i 

- 9/,+4 + 128/,+3 - 1008/,+2 + 8064/,+! 

- 14350/0/(5040 da;2), (13) 



while the mixed spatial derivatives are obtained by apply- 
ing Eq. (jlOp successively in the two direction, and have 
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the form: 

C>xy fij — 
( 9/,_4j-4 - 96/,_4j-3 + 504/,_4j-2 

- 2016/,_4j-i + 2016/,_4j+i - 504/,_4j+2 

+ 96/i_4j+3 — 9/i_4j+4 — 96/i_3j_4 

+ 1024/,_3,,-3 - 5376/,_3j_2 + 21504/,_3,,-i 

- 21504/,_3,,+i + 5376/,_3j+2 - 1024/,_3j+3 
+ 96/,_3j+4 + 504/,_2j-4 - 5376/,_2j-3 

+ 28224/,_2j-2 - 112896/,„2j-i + 112896/,_2j+i 

- 28224/,_2,,+2 + 5376/,_2j+3 - 504/,_2j+4 

- 2016/,_i,,_4 + 21504/,„ij_3 - 112896/,_ij_2 

+ 451584/,_i,,_i - 451584/,_i,,+i + 112896/,„i,,+2 

- 21504/,_i,,+3 + 2016/,_ij+4 + 2016/,+i,,^4 

- 21504/,+ij_3 + 112896/,+i,,_2 - 451584/,+i,j_i 
+ 451584/,+i,j+i - 1 12896/,+! j+2 + 21504/,+i,,+3 

- 2016/,+! j+4 - 504/,+2,j-4 + 5376/,+2,j-3 

- 28224/,+2j-2 + 112896/,+2,,-i - 112896/,+2,,+i 
+ 28224/,+2j+2 - 5376/,+2j+3 + 504/,+2,j+4 

+ 96/,+3j-4 - 1024/,+3,,-3 + 5376/,+3j-2 

- 21504/,+3j-i + 21504/,+3,,+i - 5376/,+3,j+2 

+ 1024/i-|.3j+3 — 96/i+3j+4 — 9/i+4j_4 

+ 96/,+4j-3 - 504/,+4,j-2 + 2016/,+4,j-i 

- 2016/i+4j+i + 504/j+4,j+2 - 96/i+4j+3 

+ 9 {705600 dxdy). (14) 

We modify the stencils at the refinement (and 
outer) boundary zones using the techniques proposed in 
Refs. [H^jll^]- The refinement boundary points are not 
updated during timestep, but are updated by the prolon- 
gation operation. For the first through fourth points from 
the boundary, we use standard second through eighth 
order techniques (as described below), with the excep- 
tion that we use centered derivatives for the advection 
terms if the up (down) winded derivatives do not fit on 
the grid. The first points in from the boundary are up- 
dated using standard second-order stencils, the second 
points using the standard fourth-order scheme (See [soj). 
the third points using the standard sixth-order stencils 
(See Ref. [T^), and the fourth points using the proposed 
eighth-order scheme. As in Ref. [tJ, we found satisfac- 
tory results using 6 buffer points. We use the standard 
fourth-order Kreiss-Oliger dissipation operator. 



III. RESULTS 

The initial data parameters for all new runs presented 
in this paper are given in Tables Hl HIVI For all of our 
three-black-hole runs we use a standard grid structure 
consisting of 11 levels of refinement with outer boundaries 
at 640Af and finest resolution of /i = A//80. All runs were 



TABLE L Initial data parameters for configurations with 
a third black hole intercepting a binary along the z-axis. 
(xi,yi, Zi) and {pf ,Pi) are the initial position and momen- 
tum of the puncture i, mf is the puncture mass parameter, 
mf is the horizon mass, Mf2 is the binary's orbital frequency, 
and D is the binary's initial coordinate separation. Param- 
eters not specified are zero. Configurations are denoted by 
3BHYXX, where Y = 2, 3, 4, 5 indicates the momentum of 
the third black hole, with p| = -(¥ - l)Po (See Eq. (fT^ .V 
and XX indicates the initial binary separation. 



Config 


3BH203 


3BH205 


3BH207 


yi/M 


3.6321068 


5.2079414 


6.9044853 


zi/M 


-7.2642136 


-10.415883 


-13.808971 


Pl/M 


-0.0613136 


-0.0483577 


-0.0406335 


Pl/M 
ml/M 


0.0334305 


0.0279183 


0.0242469 


0.3239234 


0.3273810 


0.3290810 


m^/M 


0.3356291 


0.3355032 


0.3351917 


y2/M 


-3.6321068 


-5.2079414 


-6.9044853 


Zi/M 


-7.2642136 


-10.415883 


-13.808971 


P%/M 


-0.0613136 


0.0483577 


0.0406335 


Pl/M 
ml/M 


0.0334305 


0.0279183 


0.0242469 


0.3239234 


0.3273810 


0.3290810 


m^ /M 


0.3356566 


0.3355216 


0.3352024 


Zi/M 


14.5284272 


20.831766 


27.617941 


Pl/M 


-0.0668610 


-0.0558366 


-0.0484938 


ml/M 


0.3247293 


0.3273813 


0.3288641 


m^ /M 


0.3313404 


0.3320131 


0.3323634 


Mfl 


0.0375000 


0.0225000 


0.0150000 


D/M 


7.2642136 


10.415883 


13.808971 



TABLE IL continuation of Table [J 

Config 3BH209 3BH211 3BH303 



performed using fourth-order techniques, except where 
otherwise noted. 



yi/M 


8.4201027 


11.117239 


3.6321068 


zi/M 


-16.840205 


-22.234477 


-7.2642136 


Pl/M 


-0.0361227 


-0.0307979 


-0.0613136 


Pl/M 
ml/M 


0.0219565 


0.0191084 


0.0668610 


0.3299485 


0.3308518 


0.3171133 


m^/M 


0.3349484 


0.3346311 


0.3299497 


y2/M 


-8.4201027 


-11.117239 


-3.6321068 


Z2/M 


-16.840205 


-22.234477 


-7.2642136 


pl/M 


0.0361227 


-0.0307979 


0.0613136 


Pl/M 
ml/M 


0.0219565 


0.0191084 


0.0668610 


0.3299485 


0.3308518 


0.3171133 


m^ /M 


0.3349593 


0.3346393 


0.3299874 


z-i/M 


33.680411 


44.468955 


14.528427 


Pl/M 
ml/M 


-0.0439130 


-0.0382167 


-0.1337220 


0.3296776 


0.3305720 


0.2955146 


m^/M 


0.3255030 


0.3327501 


0.3075078 


MQ 


0.0112500 


0.0075000 


0.0375000 


D/M 


16.840205 


22.234477 


7.2642136 
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TABLE III: continuation of Table U 



Config 


3BH307 


3BH407 


3BH507 


2/1 /A/ 


6.9044853 


6.9044853 


6.9044853 


zi/M 


-13.808971 


-13.808971 


-13.808971 


Pl/M 


-0.0406335 


-0.0406335 


-0.0406335 


Pl/M 


0.0484938 


0.0727407 


0.0969876 


ral/M 


0.3256512 


0.3197642 


0.3111199 


/M 


0.3323764 


0.3275486 


0.3206201 


2/2 /M 


-6.9044853 


-6.9044853 


-6.9044853 


Z2/M 


-13.808971 


-13.808971 


-13.808971 


Pl/M 


0.0406335 


0.0406335 


0.0406335 


P\/M 


0.0484938 


0.0727407 


0.0969876 


ml/M 


0.3256512 


0.3197642 


0.3111199 


m" /M 


0.3324005 


0.3275903 


0.3205268 


zs/M 


27.617941 


27.617941 


27.6179413 


Pl/M 


-0.0969876 


-0.1454815 


-0.1939753 


ml/M 


0.3146486 


0.2872890 


0.2319451 


mi /M 


0.3208728 


0.2987465 


0.2539505 


Mn 


0.0150000 


0.0150000 


0.0150000 


D/M 


13.808971 


13.808971 


13.808971 



A. Binary — third-black-hole interactions 

Initial data families of quasi-circular black-hole bina- 
ries in a hierarchical 3-body system with a third compan- 
ion relatively far away from the binary were studied in 
Ref. [ill . That study was based on the second-order post- 
Newtonian (2PN) approximation to the 3-body Hamilto- 
nian (We provide the 3-body Hamiltonian for our config- 
urations in Appendix [XI). In this work we use the 2-body 
Hamiltonian to generate quasi-circular black-hole binary 
configurations and then add a third black hole. This 
setup models an 'adiabatic' interaction of a binary with a 
third-black hole, where the infall timescale is significantly 
faster than the timescale for the binary to re-circularize 
through the emission of gravitational radiation. Thus 
the binary starts off with the orbital parameters it would 
have if the third body were not present. 

We evolve three sets of configurations with a third 
black hole falling towards the center of a binary perpen- 
dicular to the binary's orbital plane (which we take to be 
the xy plane). In all sets the third black hole is initially at 
a separation equal to three times the binary's separation 
and all black holes contribute equally to the total ADM 
mass given by Eq. ([7]) (As discussed below, they do not 
necessarily have similar horizon masses.). In the first set 
we choose the third black hole to have the momentum of 
a particle falling into the binary from infinity (with zero 
speed) . We obtain this momentum by assuming that the 
binary separation is fixed and then treat the resulting ef- 
fective two-body problem via Newtonian mechanics. We 
find that the third black hole has momentum 

4 / M , , 

« - ^ -^'iim- '''' 

where D is the initial binary separation and M is the 
total mass (the binary has equal momentum in the op- 



TABLE IV; Initial data parameters for the Burrau problem 
(3BHTR1 and 3BIITR2), as well as an off-center interaction 
of a binary with a third black hole falling toward the binary 
near the z-axis (3BHOC). 



Config 


3BH0C 


3BHTR1 


3BHTR2 


xi/M 


5.8677220 


5.0000000 


2.5000000 


y,/M 


4.0142579 


15.000000 


7.5000000 


zx/M 


-15.057925 


0.0000000 


0.0000000 


Pl/M 


-0.0206448 


0.0000000 


0.0000000 


pI/m 


0.0297924 


0.0000000 


0.0000000 


Pl/M 
m\/M 


0.0232357 


0.0000000 


0.0000000 


0.3298206 


0.3000000 


0.3000000 


m^ /M 


0.3355713 


0.3061460 


0.3122950 


X2/M 


-5.8750250 


-10.000000 


-5.0000000 


2/2 /M 


-4.0215595 


-5.0000000 


-2.5000000 


Z2/M 


-15.057934 


0.0000000 


0.0000000 


Pl/M 


0.0206502 


0.0000000 


0.0000000 


Pl/M 


-0.0297870 


0.0000000 


0.0000000 


Pl/M 
ml/M 


0.0232348 


0.0000000 


0.0000000 


0.3298209 


0.4000000 


0.4000000 


m" /M 


0.3355903 


0.4090649 


0.4181309 




n 0073029 


5 0000000 


2 5000000 


Vi/M 


0.0073015 


-5.0000000 


-2.5000000 


Zi/M 


30.115859 


0.0000000 


0.0000000 


Pl/M 


-5.3366e-6 


0.0000000 


0.0000000 


pI/m 


-5.3522e-6 


0.0000000 


0.0000000 


Pl/M 
ml/M 


-0.0464705 


0.0000000 


0.0000000 


0.3292338 


0.5000000 


0.5000000 


m§ /M 


0.3324495 


0.5104157 


0.5208322 


MQ, 


0.0150000 






D/M 


14.223576 







posite direction). In addition, we evolve configurations 
withp3 = ~2Po, pI = -3Fo, and pi = -4Po- We denote 
these configurations by 3BHYXX, where p| = — (Y— l)Po 
and XX indicates the binary's initial separation (XX is 
not equal to the binary separation). Initial data param- 
eters for the configurations are given in Tables Ul IIIII In 
all cases the members of the binary are given linear mo- 
menta consistent with 3PN quasi-circular orbits when ig- 
noring the contribution of the third black hole. We also 
evolve identical configurations of mass, initial positions, 
and initial momenta using a Newtonian code to com- 
pare Newtonian dynamics with the results from the fully 
nonlinear general relativistic evolutions. Note that the 
horizon mass of the third black hole is not equal to the 
horizon mass of the other two (the differences in the hori- 
zons masses between the members of the binary is a due 
solely to finite-difference errors in the isolated horizon 
algorithm [8]). This difference in horizon mass becomes 
smaller as the binary separation is increased and larger as 
the third-black-hole momentum is increased. It would be 
interesting to explore the sequence where all three black 
holes have the same horizon mass to see if qualitatively 
different results are obtained. 

In Fig. [5] we show the binary separation vector r = 
— ^2, projected onto the xy plane, for the first set of 
configurations with pg = — Poj as well as the Newtonian 
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results for these configurations (the initial third-black- 
hole to binary separation is three times the binary's ini- 
tial separation). After rescaling by the binary's initial 
separation, we see that, initially, all binaries perform sim- 
ilar (distorted) elliptical trajectories. However, the con- 
figurations with close binary separations show a prompt 
triple merger of the binary and third black hole (we indi- 
cate which systems triple merge with a * in the legend). 
As the binary's separation is increased, it undergoes more 
of an orbit before merging, and even receives a substan- 
tial kick from the third black hole. The plot seems to 
indicate that there is a critical value [i^ -Dcrit, where 
the binary quickly merges for D < Dcrit and is given a 
substantial kick above this value. The Newtonian trajec- 
tories show the opposite trend, where the larger the ini- 
tial binary separation, the more compact the orbit. The 
two systems appear to be converging to the same orbit as 
n ^ 0. The Newtonian orbits are scale invariant when 
the orbital momentum corresponds to Newtonian circu- 
lar trajectories (for the two particles in the binary). As 
is decreased and the binary separation gets large, the dif- 
ference between the post-Newtonian and Newtonian mo- 
menta becomes negligible and the system becomes scale 
invariant. To explicitly show this, we evolved similar con- 
figurations, but with MQ = IQ-^ and Mfl = 10"^", and 
confirmed that the Newtonian orbit approaches the lim- 
iting trajectory given in the figure. In Fig. [3] we show the 
position of the third black hole versus time. Note here as 
well that the Newtonian and GR trends show the oppo- 
site behaviors. In the Newtonian case, the smaller D is 
(i.e. the larger is) the more likely the third particle is 
to pass through the binary, while in the GR simulations 
the third black hole always merges with the binary for 
small D. The GR and Newtonian evolutions approach 
each other as D tends toward infinity. 

In Figs. [4] and [5] we explore how the three-body sys- 
tem behaves when the initial z-momentum is increased, 
while still keeping the initial third-black-hole to binary 
separation at three times the binary's initial separation. 
We choose a relatively large initial binary separation of 
D = 13.808971M = 20. 7134565Mb (Mb is the binary's 
mass). In the absence of the third black hole, this binary 
would complete approximately 25 orbits before merging. 
In Newtonian theory, configurations 3BH207, 3BII407, 
and 3BH507 (i.e. = -Po, -3Po, -4Po) result in a 
bound binary with the third particle ejected to infinity, 
while for configuration 3BH307 the entire system is dis- 
rupted; leading to three free particles. Note also that 
in Newtonian theory the third black hole passes through 
the binary and is ejected toward z = —oo for all configu- 
ration except 3BII207. The GR simulations on the other 
hand seem to indicate that the binary is disrupted for all 
configuration except 2BII207 (which triple merges). In- 
terestingly the behavior of the third particle approaches 
the Newtonian behavior as is increased (as can be seen 
in Fig. ED. Note that, in the GR simulation of 3BH307, 
the third black hole was not able to escape to z = — oo, 
but rather was bounced toward z = -|-oo. 



FIG. 2: The xy projection of the binary-separation trajectory 
when the third black hole falls towards the center of the binary 
along the z axis with initial momentum pi = — Po- The co- 
ordinates have been rescaled by the initial binary separation. 
The solid curves are the GR trajectories, while the dot-dashed 
curves are the Newtonian trajectories. Mimit is the Newto- 
nian trajectory for D —f oo (i.e. ^ 0). There seems to be 
a critical separation between 3BH207 and 3BH209 where the 
system transitions from a prompt triple merger to an ellipti- 
cal binary plus third black hole (we indicate which systems 
triple merge with a * in the legend). 
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We conclude our analysis of this three-black-hole con- 
figuration by comparing the behavior of the system when 
the z-momentum is doubled and the binary separation is 
increased (while still keeping the initial third-black-hole 
to binary separation at three times the initial binary sep- 
aration). In Fig. [5] we show the xy trajectories for config- 
urations 3BH203, 3BH303, 3BH207, 3BH307. Note that 
for the closer binary (3BH203, 3BH303) increasing the 
momentum of the third black hole (which increases the 
eccentricity of the binary) causes the binary to merge 
sooner, while for the larger binary this same increase 
leads to a (possible) disruption of the binary. 

It is interesting to compare the dynamics of this three- 
black-hole system under modified initial conditions that 
reduce the symmetry of the problem. We began this 
type of analysis by considering a purely Newtonian sys- 
tem consisting of a circular binary with orbital frequency 
Mn = 0.015 and mass Mb = 2/3M (similar to 3BHY07, 
but with Newtonian, rather than post- Newtonian, orbital 
momenta) and a third particle of mass 1/3M located 
at {x,y,z)/D = (0.0099,0.0099,1000). We then evolved 
this configuration using purely Newtonian evolution un- 
til the third-particle — binary separation was ~ 3-D. We 
then took the Newtonian position and momentum pa- 
rameters and evolved them using full GR. In the purely 
Newtonian evolution, this 3-body system leads to an ex- 
change of partners in the binary, where particle 2 is 
ejected from the system and particles 1 and 3 form a 
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FIG. 3: The z-component (rescaled by the initial binary sep- 
aration) of the trajectory of the third black hole with initial 
momentum p| = — Pq. Q is the binary's initial frequency. 
The solid curves are the GR trajectories, while the dot-dashed 
curves are the Newtonian trajectories. A^umit is the Newtonian 
trajectory for D —f oo (i.e. f2 — > 0). Systems that triple-merge 
are indicated by a * in the legend. 

4 




FIG. 5: The z-component (rescaled by the binary initial 
separation) of the trajectory of the third black hole with initial 
momentum = —Po, —2Po, — 3Po, — 4Po. fl is the binary's 
initial frequency. The solid curves are the GR trajectories, 
while the dot-dashed curves are the Newtonian trajectories. 
Systems that triple-merge are indicated by a * in the legend. 
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FIG. 4: The xy projection of the binary-separation tra- 
jectory when the third black hole falls towards the cen- 
ter of the binary along the z axis with initial momentum 
P3 = — Po, — 2Po, — 3Po, — 4Po. The coordinates have been 
rescaled by the initial binary separation. The solid curves 
are the GR trajectories, while the dot-dashed curves are the 
Newtonian trajectories. Note that 3BH207 and 3BH507 re- 
sult in bound binaries in Newtonian theory. There appears to 
be a critical momentum where the system transitions from a 
prompt merger to a highly-elliptical, or perhaps even hyper- 
bolic, orbit. Systems that triple-merge are indicated by a * 
in the legend. 
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FIG. 6: The xy projection of the binary-separation trajectory 
when the third black hole falls towards the center of the bi- 
nary along the z axis with initial momentum = — Po, — 2Po 
at two different initial binary separations. The coordinates 
have been rescaled by the initial binary separation. Here the 
critical separation between prompt-merger and (possible) dis- 
ruption is a function of both the binary separation and the 
initial z-momentum. Systems that triple-merge are indicated 
by a * in the legend. 
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new binary with eccentricity e = 0.645132 (See Fig. [7]). 
The full GR simulations, on the other hand, show a triple 
merger that merges faster than 3BH207 (which likely due 
to our using Newtonian, rather than post-Newtonian, or- 
bital momenta). In Fig. [5] we show the xy projection 
of the initial binary trajectory for 3BH0C, as well as 
3BH207 and the Newtonian trajectory. 



8 



FIG. 7: The xy projection of the Newtonian trajectories for 
3BHOC. The initial positions of the three particles are indi- 
cated by filled circles. Note that the initial P1-P2 binary is 
disrupted and a new P1-P3 binary is formed (The new P1-P3 
binary is inclined with respect to the xy plane.). The initial 
motion of P3 is essentially along the 2-axis and is thus not 
apparent. 
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FIG. 9: The fuU Newtonian tracks for the 3BHTR1 configu- 
rations (3BHTR2 is obtained by rescaling by a factor of 1/2). 
Note that particle 1 is ejected and moves towards the upper 
right of the figure, while particles 2 and 3 form an elliptical 
binary that moves towards the lower left. Initial positions are 
indicated by a filled circles. Particle 1 (mass 0.3M) is the one 
on top, particle 2 (mass QAM) is the lower left, and particle 
3 (mass 0.5M) is the lower right. The arrows indicate the 
trajectories of the recoiled binary and lone particle. 




FIG. 8: The xy projection of the binary-separation trajectory 
for 3BHOC, as well as 3BH207 and the Newtonian trajectory. 
Note that 3BH0C is more elliptical (and thus merges sooner) 
than 3BH207. 
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B. Burrau Three-Body Configuration 



FIG. 10: The GR and Newtonian tracks for configurations 
3BHTR1 and 3BHTR2, the latter rescaled by a factor of 2 (the 
Newtonian trajectories for these configurations are identical 
up to a scaling). After rescaling, the two GR trajectories 
are very similar (see right inset), and differ significantly from 
the Newtonian trajectories after the first BH2 — BH3 close 
encounter (see left inset), which leads to a merger in GR. 
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Configurations 3BHTR1 and 3BHTR2 belong to a set 
of configurations known as the Burrau three-body prob- 
lem J79,, iS^, This system consists of three particle, 
initially at rest, arranged at the vertices of a right trian- 
gle with sides of length 3p, 4p, 5p and masses 4/.t, 5/i, 
where the particle of mass i ^ is located on the vertex op- 
posite the side of length i p. After suitable rescaling, the 
Newtonian trajectories are independent of p and [i. In 
Newtonian theory, this configuration will lead to particles 
2 and 3 forming a highly elliptical binary (with elliptic- 
ity e ~ 0.989), and particle 1 and the 2-3 binary ejected 



in opposite directions (See Fig. [H]). The trajectories of 
the corresponding black holes in a full GR simulation are 
quite different (although, if p is kept fixed and p is in- 
creases, eventually the GR simulations will reproduce the 
Newtonian trajectories). Here we see that the trajecto- 
ries scale reasonably well with p (i.e. compare 3BHTR1 
and 3BHTR2 in Fig. [rU]) . but rather than forming a bi- 
nary and free particle, this system quickly merges to form 
a single black hole. 
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FIG. 11: The = 2,m = 2) mode of V'4 for 3BHTR1 showing 
both the two mergers. Note that the imaginary part of this 
mode is more sensitive to the second merger. The two mergers 
are clearly seen in the absolute value of '04 • 
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FIG. 12: The = 2, m = 2) mode of ^4 for 3BHTR2. Here 
the real part of i/)4 seems to only indicate the presence of one 
merger, while the second merger is more clearly visible in the 
imaginary part. 
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The Newtonian, and even 2PN, trajectories for this 
three-black-hole problem will not agree with the GR tra- 
jectories if the closest approach of any two BHs is within 
about ten times the combined mass of the two BHs. To 
analyze the behavior of the system at different scales, we 
set the masses of the three holes to 3m, Am, 5m, and set 
the initial separations to ri3 = 20Qmp, r23 = 150mp, 
ri2 = 250mp. The initial merger of BH2 and BHS 
seen in the GR simulation corresponds to a close ap- 
proach of r23/(TO2 + "^3) = 0.054p. Thus, for the ini- 
tial close encounter not to lead a quick merger, we need 
p ^ 185, which means that the initial value of will 
need to be larger than r23 = 1.4 x 10~Vc("i/Mq). How- 
ever, the closest approach over the entire trajectory is 
^12/ {rrii + '713) = 0.0023p, and in order for this approach 
not to lead to a prompt merger the initial value of r23 
would need to be r23 = 3.3 x lQ~^pc{m/MQ). 



FIG. 13: The = 2, m = 2) mode of i/'4 for Q38 using both 
4th and 8th-order algorithms. Note that the phase error in 
the M/80 8th-order waveform is apparently smaller than the 
phase error in the M/120 4th-order waveform. 
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C. Eighth-Order Accuracy 



The waveforms from 3BHTR1 and 3BHTR2 show dou- 
ble mergers, with the two mergers more closely spaced in 
3BHTR2. In Figs. [Hand [H we show the = 2,to = 2) 
mode of ?/'4 extracted at r = 40M = 33.33Madm (here 
Madm — 1.2M). Interestingly, the real part of the 
{£ — 2, m = 2) mode is more sensitive to the initial 
merger (of BHs 2 and 3), while the imaginary part is 
more sensitive to the second merger (of the 2-3 remnant 
with BH 1). The radiated energy and angular mo- 
mentum for configuration 3BHTR1 were i?rad/.^^ADM = 
(8.5 ± 1.0) X 10-'' and Jrad/A/ADM = (2-8 ± 2.1) x lO"'', 
while for configuration 3BHTR2 the radiated quantities 
were Er^^/M adm = (6.6 ±0.2) x lO"* and J^ad/M^DM = 
(1.2 ±0.6) X 10"''. The recoil velocities for 3BHTR1 and 
3BHTR2 were (4.1 ± 2.0)kms~i and (3.0 ± 0.6)kms-i 
respectively. 



We evolved several configurations, including Q38 of 
Refs. [13, [13] (which is a quasi-circular, non-spinning 
binary with mass ratio 3:8) and 3BH102 of Ref. [53| 
(which is a three-black-hole configuration with planar 
orbits), using our standard fourth-order code and the 
new eighth-order code. In all cases we used the same 
grid structure with a maximum resolution of M/80. In 
Fig. [T3] we plot the (i = 2, m = 2) component of V'4 for 
the Q38 configuration with a Gamma-driver parameter 
of ?7 = 2/AI. As pointed out in ^33\, this choice of rj leads 
to a low effective resolution for our lowest resolution run 
(i.e. M/80), leading to a large phase error. The corre- 
sponding eighth-order run appears to be more accurate 
than the A//120 fourth-order run (i.e. its phase is closer 
to that predicted by extrapolating using the fourth-order 
Af/100 and A//120 runs). Eighth-order also significantly 
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FIG. 14: The z component of the trajectory of BHl for con- 
figuration 3BH102. This component should be zero by sym- 
metry. The 8th-order algorithm produces an error that is 3.33 
times smaller. Note that the error in the 4th-order trajectory 
is smaller than the central resolution of M/80 until t ~ 290. 
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FIG. 15: The xy projection of the trajectory of all three BHs 
in 3BH102. Solid curves were produced using the fourth- 
order algorithm; dot-dashed curved using the eighth-order al- 
gorithm. 
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reduces the error in three-black-hole simulations. To test 
this, we evolved configuration 3BH102 (which consists 
of three equal-mass black holes, labeled BHl - BHS in 
Fig. [ini in a planar orbit) using a shifted grid structure 
(i.e. one not centered on the origin) that does not exhibit 
the z-reflection symmetry of this configuration. In Fig. [HI 
we plot the z-component of the trajectory for BHl (which 
should be identically zero). Note that late-time error in z 
is reduced by a factor of 3.33 when going from 4th to 8th- 
order for this coarse grid structure, which may be crucial 
for longer term evolutions. The xy trajectories do not 
change qualitatively for 3BH102 when moving to eighth- 
order, as is apparent in Fig. [151 We will present results 
from convergence and efficiency tests in a forthcoming 
paper. Here we note that we found fourth-order conver- 
gence (See Fig. [TB)) in the waveform for three-black-hole 



FIG. 16: The real part of the__{£ = 2, m = 2) mode of ^4 
for configuration 3BH1 of Ref. 58] with central resolutions of 
Af/80, M/96, and A'//115.2, along with a convergence plot of 
data. Note the fourth-order convergence, and the smallness 
of the waveform amplitude. 
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systems using our standard algorithm. 



IV. DISCUSSION 

Our fully nonlinear evolutions of three-black-hole sys- 
tems demonstrate dissipative General Relativistic effects 
due to emission of gravitational radiation and black- 
hole mergers, which do not have Newtonian (or even 
2PN) counterparts, that dramatically change the qual- 
itative behavior of the system. These effects become 
important when two objects approach each other closer 
than about 100 times their combined mass [s^ (or d ~ 
100G(mi -f ra2)l<? in more conventional units) and be- 
come dominant when the two objects approach within 
about 10 times their combined mass. Even when the ef- 
fect is small, the sensitivity of three (and more) body en- 
counters to small perturbations makes including PN and 
GR effects important for obtaining the correct dynamics 
of an N-body simulation. The GR effects in particular 
provide a natural regularization of the Newtonian prob- 
lem on small scales. 

An important problem in galactic dynamics is deter- 
mining how the merger process of the two (initially) cen- 
tral black holes proceeds after the collision of the host 
galaxies [s^ . [ssj . In our simulations we confirm that 
the resulting binary will, in general, become highly ellip- 
tical if it interacts with a similar mass black hole. These 
eccentric orbits will drive the binary merger due to emis- 
sion of gravitational radiation. However, we also found 
that the interaction can produce an immediate merger of 
the triple system. If we set physical scale of our simula- 
tion by setting the black-hole masses to lO^M©, then the 
prompt mergers that we observe occur at separations of 
the order of milli-parsecs. These three-body interactions 
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provides a possible mechanism for resolving the "last par- 
sec problem." [HI We will present results from the ef- 
fects of unequal-masses, spins, and further separated bi- 
naries in a forthcoming paper, where we will also explore 
the potential for critical phenomena. It is important to 
note that these three-black-hole interactions provide a 
mechanism for producing highly-elliptical close-binaries 
(which would otherwise have circularized due to emis- 
sion of gravitational radiation during the inspiral) like 
those studied in Ref . 1431 . 
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APPENDIX A: 2PN ORBITS 

In this paper we did not include all Newtonian or Post- 
Newtonian 3-body interactions when producing initial 



data for quasi-circular orbits. However, these interaction 
can be taken into account, up to 2PN, using the Hamil- 
tonian of the 3-body interactions in the ADM gauge as 
given by Eq. (5) in Ref. [13 (and Eqs. (l),(2),(Al) in 
Ref. H). 

Our configurations are characterized by 



mi 

Xi = 

X2 = 

2:3 = 

pI-~ 



- ra-i = 7713 = m 
0,7/1 = r/2,zi = -r 
0,7/2 = -rl2, Z2 = -r 
0, 7/3 = 0, Z3 2r 
-J/r,p\^0,pl^-pl/2 
J/r,pl^Q,pl = -pl/2 
0,p| = 0, 




(Ala) 
(Alb) 
(Ale) 
(Aid) 
(Ale) 
(Alf) 
(Alg) 

(Alh) 



where r is the separation of the binary, the initial third- 
body — binary separation is 3r, and the initial linear mo- 
mentum of the third body corresponds to the approxi- 
mate momentum it would have after falling toward the 
binary from infinity. 

The process for computing the quasi-circular orbits 
proceeds as follows. We first choose a value of the orbital 
angular momentum, J, and then compute the linear mo- 
menta of the holes as given above in Eqs. (|A1[) . We then 
compute the Hamiltonian, 



H = (3241792 (l9 + V2^\ + 21904r (4 (-1369 + 269\/37) 777^ - 37 (l9 + VW) (87^7^ - 3pf) r) + 



148r2 (1152 (90781 -I- 8503\/37) 777^-8 (296 (9583 + 709^37^ 777^ -|- (l3801 -I- 19075^37^ pf ) rra^ + 
1369 (19 + A/37) (128777" - 16pf 777^ + 3pf ) r^) + (-128 (9888649 + 1247443^37^ 



Q 

777 - 



9472 (74 (1079 4- I25V37) 7^7^ + 3 (100109 + 6I40V37) pf ) rm^ - 4 (l75232 (851 + 113%/37) 777"-^ 
592 (143449 + 85651\/37) pf 777^ + 3 (26603 - 429895\/37) pf ) r^m^+ 

151959 (19 + V37) pf (1287^7" - 24pf 777^ + llpf ) r^) )/(25934336 (l9 + VS?) 777^^), (A2) 



and find the value of r that gives a local minimum (i.e. the binary in a quasi-circular orbit. Some of the results 
drH ~ 0) [ill, which gives the value of the separation of of this process are summarized in Table IVl 
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